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<— 1 1 ABSTRACT 

> : 

Context. Massive young stellar objects (MYSO) are surrounded by massive dusty envelopes, whose physical structure and geometry are 
determined by the star formation process. 

Aims. Our principal aim is to establish the density structure of MYSO envelopes on scales of ~ 1000 AU. This constitutes an increase of a factor 
~10 in angular resolution compared to similar studies performed in the (sub)mm. 

Methods. We have obtained diffraction-limited (0.6") 24.5 fim images (field of view of 40" x 30") of 14 well-known massive star formation 
regions with the COMICS instrument mounted on the 8.2 meter Subaru telescope. We construct azimufhally averaged intensity profiles of the 
resolved MYSO envelopes and build spectral energy distributions (SEDs) from archival data and the COMICS 24.5 pm flux density. The SEDs 
range from near-infrared to millimeter wavelengths. Self-consistent 1-D radiative transfer models described by a density dependence of the 
^» , form n{r) cc r~ p are used to simultaneously compare the intensity profiles and SEDs to model predictions. 
J^j ' Results. The images reveal the presence of discrete MYSO sources which are resolved on arcsecond scales, and, to first-order, the observed 
rS , emission is circular on the sky. For many sources, the spherical models are capable of satisfactorily reproducing the 24.5 pm intensity profile, 

■ the 24.5 yum flux density, the 9.7 pm silicate absorption feature, and the submm emission. They are described by density distributions with 
p = 1.0 ± 0.25. Such distributions are shallower than those found on larger scales probed with single-dish (sub)mm studies. Other sources have 
density laws that are shallower/steeper than p = 1.0 and there is evidence that these are viewed near edge-on or near face-on respectively. In 
these cases spherical models fail to provide good fits to the data. The images also reveal a diffuse component tracing somewhat larger scale 
structures, particularly visible in the regions S 140, AFGL2136, IRAS 20126+4104, Mon R2, and Cep A. 

Conclusions. We find a flattening of the MYSO density law going from scales probed with single-dish submm observations down to scales 
of ~1000 AU probed with the observations presented here. We propose that this may be evidence of rotational support of the envelope. This 
finding will be explored further in a future paper using 2-D axisymmetric radiative transfer models. 
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1. Introduction 

The early stages in the lives of massive stars are identi- 
fied with highly luminous yet very cool objects deeply em- 
bedded in molecular clouds. Their emission is characterised 
by a steeply rising spectral energy distribution (SED) peak- 
Send offprint requests to: W.J. de Wit, e-mail: 
w. j .m.dewit@leeds.ac.uk 

* Based on data collected at the Subaru telescope, which is operated 
by the National Astronomical Observatory of Japan. 



ing around 100 ^m, which features strong silicate absorption. 
These properties are indicative of radiation reprocessed by a 
dusty circumstellar envelope. Measured bolometric luminosi- 
ties are such that if a single main sequence star resides at the 
heart of the dusty envelope, it should have the ability to ionise 
its surroundings, yet only little (if any) recombination radia- 
tion is observed. This radio-quiet appearance of massive young 
stellar objects (MYSO) is unlike that of ultra-compact (UC) 
H n regions (e.g. Hoare et al. 2007), and the latter can therefore 
be considered as a successor phase. The reason may be that 
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the MYSO is actively accreting material from its surrounding 
environment quenching the development of an H n region. The 
ongoing accretion would also give rise to bipolar outflow activ- 
ity ubiquitously observed in massive star forming regions. It is 
clear that MYSOs are prime observational targets for address- 
ing outstanding issues in our current understanding of massive 
star formation. 

In this paper we address the radial structure of the MYSO 
dust envelope. It is determined by the forces that operate during 
the onset and the subsequent evolution of the initial molecular 
core. The radial density profile is predicted to be a powerlaw 
with a value for the power index that depends on the domi- 
nant physics. The exact power index can be extracted from the 
observables by the use of radiative transfer models. Spherical 
envelope models may be assumed for the dust that dominates 
the emission at wavelengths larger than 30 yum. At these wave- 
lengths the SED of MYSOs (but also UCH us) are observed to 
be remarkably similar, arguing for little deviation from spheri- 
cal symmetry (Chini et al 1986; Churchwell, Wolfire & Wood 
1990;Henningetal 1990; Hoare et al. 1991; Giirtleret al. 1991; 
Wolfire & Churchwell 1994; Faison et al. 1998; Hatchell et al. 
2000; Van der Tak et al. 2000; Mueller et al. 2002; Beuther et 
al. 2002). At shorter wavelengths (A < 30//m) the geometry 
of the envelope becomes important (Yorke & Sonnhalter 2002; 
Whitney et al. 2003). For example, under favourable inclina- 
tions, mid-IR radiation can be observed to originate directly 
from the surface of cavity walls that are sculpted by the polar 
outflows. In this case, the mid-IR photons are emitted by warm 
dust particles that have a clear line-of-sight to the star (e.g. De 
Buizer 2007). At even shorter wavelengths, near-IR photons 
from the (generally monopolar) reflection nebulae (e.g. Alvarez 
et al. 2004, 2005) may originate either from the stellar surface, 
an inner dust truncation structure or from an accretion disk. 
They can scatter and escape through existing inhomogeneities 
in the spherical envelope (e.g. Giirtler et al. 1991; Henning et 
al. 2000) and still suffer extinction from any foreground molec- 
ular cloud material (e.g. De Buizer, Osorio & Calvet 2005). 

Here, we aim to constrain the radial density distribution 
on scales of 1000AU using resolved 24.5 /im emission. This 
constitutes an increase of a factor 10 in angular resolution 
compared to similar studies performed in the (sub)mm. We 
present diffraction-limited images at 24.5 //m, the longest mid- 
IR wavelength amenable to high resolution imaging from the 
ground with large telescopes. This long mid-IR wavelength 
maximises the possibility of resolving the envelope emission 
because, due to the nature of the temperature gradients in the 
optically thick emitting region, the size of the emission re- 
gion gets larger with increasing wavelength to the power of 
about 1.5, whilst the diffraction limit of the telescope only in- 
creases linearly with wavelength. We have selected a set of 14 
well-known MYSO and imaged these at 24.5 /urn with the 8m 
Subaru telescope. The images have an angular resolution set 
by the telescope's diffraction limit of 0.6", corresponding to 
linear scales of ~1000AU for the average distance of 1.5 kpc 
to our target MYSOs. Most previous 10/Lmi and 20/j.m imag- 
ing of MYSOs have been carried out on 4 m class telescopes 
where the (radio-quiet) objects are invariably unresolved (e.g. 
Mottram et al. 2007; De Buizer et al. 2000, 2005). Some work 



Table 1. COMICS 24.5 /im observing log of standard stars (P = PSF 
standard; S = standard). Standard stars are further divided into Cohen 
(C) and Sloan (L) standards. Upper part of the table are observations 
performed with the Q24.5-OLD filter, the lower part denote the ob- 
servations performed with the Q24.5-NEW filter. Filter characteristics 
are given in Table 3. 
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Asteroid #5 1 1 (P) 
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1.8 
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2002/12/15 


Asteroid #51 (P) 
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2.0 
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2003/06/17 


fi Cep (SL) 


402 


1.3 


1800 


2003/06/20 


li Cep (SL) 


402 


1.3 


1000 


2003/07/14 


fi Cep (SL) 


101 


1.5 


180 


2003/11/12 


a Tau (SL) 


1204 


1.2 


800 


2004/01/08 


a Tau (SL) 


402 


1.0 


830 


2004/05/05 


H Cep (SL) 


201 


1.4 


790 


2004/06/08 


a Her (SL) 


100 


1.2 


630 


2004/06/08 


H Cep (SL) 


602 


1.3 


3000 


2004/06/08 


a Sco (SL) 


201 


1.9 


1200 


2004/07/14 


H Cep (SL) 


401 


1.3 


1500 


2005/07/27 


a Her (SL) 


201 


2.6 


160 


2005/07/27 


yAql(SC) 


201 


1.3 


30 


2005/12/13 


a Tau (SL) 


401 


1.3 


430 


2005/12/21 


a Tau (SL) 


1203 


1.4 


620 



on 8-10m class telescopes at 8-22 //m has begun to resolve a 
few sources such as BN and source n in Orion (Shuping et al. 
2004; Greenhill et al. 2004), but no detailed modelling of these 
data has been carried out. 

We analyse the resolved emission in conjunction with the 
SED in terms of spherical dust radiative transfer models as 
calculated by DUSTY, and using background literature in- 
formation for each individual case. We present simultaneous 
model fits to the 24.5 /im intensity profile and the SED, that 
stretches from the near-IR to the mm wavelengths. It includes 
the silicate 9.7/im absorption profile thanks to ISO-SWS spec- 
tra for nearly the whole sample. The approach of simultane- 
ously analysing the intensity profiles and SEDs follows van der 
Tak et al. (2000), Hatchell et al. (2000), Beuther et al. (2002), 
Mueller et al. (2002), Hatchell & van der Tak (2003), Williams 
et al. (2005). Most of the quoted work is exclusively aimed at 
(sub)mm wavelengths, and thus probing linear scales about a 
factor 10 larger than in the mid-IR. Constraints on the density 
distributions from (sub)mm and mid-IR are therefore highly 
complementary as they give insight into the evolution of the 
density distribution as function of radius. 

This paper is organized as follows. Our observational data 
were taken with the Japanese 8.2 meter Subaru telescope on 
Hawaii in conjunction with the COMICS instrument. We give 
an overview of the instrument and detail the observations in 
Sect. 2. We present and discuss the morphology seen on the 
COMICS images in the subsections to Sect. 3.3 and Sect. 3.4, 
alongside the simultaneous modelling of the intensity profile 
and SED. We summarise the modelling part highlighting cer- 
tain trends in Sect. 4. We discuss our findings and their conse- 
quences in Sect. 5 and conclude in Sect. 6. 
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Table 2. COMICS 24.5 /im observing details. Integration is the duration of on-source integration in seconds, and AM is the average airmass. The 
signal to noise in Col. 5 is calculated as peak flux divided by the background standard deviation. Observations after 2004/06/08 are performed 
with Q24.5-NEW filter. The J2000 coordinates in Col. 6 and 7 correspond to the target MYSO in each image, with its reference given in Col. 8. 
The offsets in Col. 9 correspond to the identification of known sources (Col. 10) in the images. The one uncertain identification is marked with 
a colon. The uncertainty of the flux densities given in Col. 1 1 is on the order of 10%. The Mon R2 and W3 regions are relatively large and have 
been mosaiced using five, respectively four images (two long and two short exposures). 
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(0,18) 
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M8E 
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18:04:53.3 
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(0,0) 

(-6,5) 
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IR 
HII 

MIRSl 


210 
30 
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AFGL2136 


2005/07/27 


201 


2.5 


170 


18:22:26.5 


-13:30:12.0 


3 


(0,0) 


IR 


140 


AFGL2591 


2004/05/05 


603 


1.1 


2000 


20:29:24.9 


+40:11:20.3 


4 


(0,0) 

(-5,-3) 

(-7,9) 


IR 
HII 

MIRSl 


870 
170 
20 


NGC2264 


2002/12/15 


1306 


1.9 


1100 


06:41:10.1 


+09:29:34.0 


5 


(0,0) 


IRSl 


330 


S255 


2003/11/12 


1404 


1.0 


1000 


06:12:54.1 


+ 17:59:25.1 


6 


(0,0) 
(-2.5,0.5) 


IRS3 
IRSl 


170 

20 


AFGL5180 


2005/12/20 


1003 


1.2 


100 


06:08:53.3 


+21:38:30.5 


4 


(0,0) 
(12,-4) 

(2,2) 


IRSl 
HII 

MIRS3 


490 
210 

35 


IRAS 20 126 


2005/07/27 


401 


1.3 


60 


20:14:26.1 


+41:13:32.5 


7 


(0,0) 


IR 


60 


Mon R2 


2005/12/13 


100 


1.2 




06:07:47.8 


-06:22:54.7 


1 


(0,0) 

(-31,3) 

(-33,17) 


IRS3 
IRS2 
IRS5 


1150 
40 
70 


AFGL437 


2005/12/21 


401 


1.3 


80 


03:07:24.6 


+58:30:44.4 
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(0,0) 

(0,10) 

(-6,6) 


S 
N 
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30 
30 
200 


AFGL4029 


2005/12/21 


1203 


1.3 


80 


03:01:31.3 


+60:29:12.9 
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(0,0) 
(24,1) 


IRSl 
IRS2 


70 
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AFGL961 


2003/11/12 


1204 


1.0 


620 


06:34:37.7 


+04:12:44.4 


10 


(0,0) 
(-5,-2) 


E 
W 


250 
60 
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2005/12/21 


501, 100 


1.4 




02:25:41.4 


+62:06:21.8 
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(0,0) 
(3,-7) 
(-1,-12) 
(-2,11) 


IRS5 
IRS6 
IRS7 
MIRSl 


1300 
70 
160 
15 


Cep A 


2004/07/13 


602 


1.4 


70 


22:56:18.0 


+62:01:49.5 


12 


(0,0) 




440 



(1) Hackwell et al. (1982); (2) Simon et al. (1984); (3) Kastner et al. (1992); (4) Tamura et al. (1991); (5) Thompson et al. (1998); (6) 
Longmore et al. (2006); (7) Sridharan et al. (2005); (8) Wynn- Williams et al. (1981); (9) Zapata et al. (2001); (10) Castelaz et al. (1985); (11) 
van der Tak (2005); (12) comparison with Spitzer images, see Fig. 29 



2. Observations and data reduction 

2.1. Observations with the COMICS instrument 

Tables 1 and 2 summarise the 24.5/im observations of the 14 
target MYSOs and the standard stars. All measurements were 
made using the mid-infrared imaging spectrometer COMICS 
(Kataza et al. 2000) at the Cassegrain focus of the 8.2 meter 
Subaru Telescope on Mauna Kea, Hawaii. Imaging mode of 
COMICS utilises a Raytheon 320x240 Si:As IBC array, which 
is cooled by a Sumitomo 4-K Gifford-McMahon type cryo- 
cooler but usually operates at around 7~8 K because of the self- 
heating. The camera provides over-sampled diffraction-limited 



images at 24.5 //m with a pixel size of 0.13x0.13 arcsec 2 and a 
field of view of approximately 40x30 arcsec 2 . 

We used two different 24.5 fim filters which are both man- 
ufactured by Infrared Multilayer Laboratory, University of 
Reading. The characteristics of these filters are summarised 
in Table 3, and the transmission curves, along with an atmo- 
spheric transmission model above Mauna Kea, are shown in 
Figure 1. As can be clearly seen in the plot, the new filter 
(Q24.5-NEW) is a much better fit to the small atmospheric win- 
dow at 24.5 yum. As a result, the whole array can be read out 
with the new filter, whilst only a part of the array is read out 
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Fig. 1. Transmission curves for the Q24.5-OLD filter (solid line), 
Q24.5-NEW (dashed line), and the atmosphere above Mauna Kea 
(dotted line). The filter transmission spectra have been made available 
by the manufacturer. The atmospheric transmission curve has been 
calculated using USF HITRAN-PC for a standard 'US tropical model' 
with an H2O partial pressure of 1.35 mmHg at 4200 m looking through 
the atmosphere at a zenith angle of 30° (Tomono 2000). 



with the old one since radiation from the sky quickly saturates 
the well. 

Chop-subtracted frames were stacked using the imalign 
task (with cubic-spline interpolation) in the iraf 1 data reduction 
package. Flux calibration was achieved against either Cohen 
(Cohen et al. 1995; 1999) or Sloan (Sloan et al. 2003) stan- 
dards. We selected stars at similar airmasses to the target 
MYSOs whenever possible; however, in some cases when this 
was not feasible, we scaled the standard flux to the airmass 
of the relevant object by the atmospheric extinction relation- 
ships measured on 2002/12/15 UT (0.57 mag per airmass) and 
on 2003/11/12 UT (0.56) for the Q24.5-OLD filter. Note that 
these extinction values would probably only apply to this spe- 
cific filter, along particular lines of sight, and at the times of ob- 
servations. In some other instances (e.g. 2005/12/20 UT), even 
this airmass correction was not possible due to lack of appro- 
priate data, and we reluctantly accepted the airmass mismatch 
as additional uncertainty. We estimate the overall uncertainty in 
flux calibration to be of the order of 10%. The largest contribu- 
tion usually comes from the absolute calibration uncertainty in 
the standard flux templates. Details of the target observations 
and extracted fluxes are given in Table 2. Images are not astro- 
metrically calibrated. Positional offsets of the various sources 
in our images are with respect to the brightest 24.5 yum source, 
generally identifiable with the brightest MYSO in each region 
presented in this paper. 

2.2. Point spread function reference stars 

We spend a few words on the point spread function (PSF) 
reference objects, as they are crucial in our analysis of com- 



1 iraf is distributed by the National Optical Astronomy 
Observatories, which are operated by the Association of Universities 
for Research in Astronomy, Inc., under cooperative agreement with 
the National Science Foundation. 
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Fig. 2. Comparison of azimuthally averaged intensity profiles of the 
observed standard stars. Standard stars represented by open symbols 
are apparently surrounded by extended envelopes (see text). 

Table 3. Characteristics of two 24.5//m filters. AA is measured from 
the 50% transmission cut-on to cut-off wavelengths and A is the half- 
way point between the two. Measurements made available by the man- 
ufacturer. 



ID 


^0 


AA 


Peak transmission 




Oum) 


w 


(%) 


Q24.5-OLD 


24.47 


1.91 


64 


Q24.5-NEW 


24.56 


0.75 


64 



paring model images to the resolved MYSO emission (see 
Sect. 3.3). All standard star observations are given in Table 1. 
To test the validity of both the PSF standards and the flux 
standards as a model for the instrumental PSF, we investigate 
their azimuthally averaged intensity profiles in Fig. 2. The fig- 
ure makes clear that three of the flux standard stars are actually 
extended, and therefore do not qualify as a PSF standard. They 
are the following objects: 

// Cep is one of the largest stars in our Galactic neighbourhood. 
The star was observed a total of six times with COMICS. 
Detailed analysis shows this star to be extended in such 
degree and complexity, that we report on this star in a separate 
publication (de Wit et al. 2008). 

a Sco is a supergiant star undergoing mass-loss. Mid-IR 
images at 12.5/im and at 20.8/im clearly show an extended 
distribution of circumstellar dust (Marsh et al. 2001). 
a Her is an irregular variable M star. Although the star 
mimics a PSF in the inner part of the profile (Fig. 2), its extent 
becomes apparent at radii > 2". Interestingly, interferometric 
observations of the star have been interpreted as showing a 
mass loss event in the period 1989-1992. The material has an 
expansion speed of approximately 170 mas yr~' (Tatebe et al. 
2007), and it should have reached a distance of 2.5" from the 
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star at the time of the observations. This distance is comparable 
to the distance where the excess emission becomes apparent. 

As a consequence of the above discoveries, most observa- 
tions do not have a concomitant PSF reference observation. In 
addition, the instrumental PSF shows a non-symmetric pattern 
due to the mirror support structure that rotates on the sky, mak- 
ing the interpretation of non-symmetric structure in the par- 
tially resolved sources taken at different times difficult. A com- 
parison of the azimuthally averaged intensity profile of the six 
fi Cep observations taken over a period of more than 1 year 
shows however very few differences (see de Wit et al. 2008). 
This provides confidence in the stability of the (azimuthally 
averaged) COMICS PSF, and we therefore construct a single 
reference PSF from the remaining three genuine point sources 
(a Tau, and asteroids 5 1 and 511), which we use throughout the 
remainder of this paper. 

3. Results 

3.1. Brief description of the images 

The 24.5 pm images reveal that all principal MYSO sources in 
the fourteen targeted massive star forming regions are resolved. 
The sources are generally discrete and have fairly simple, cir- 
cular morphologies on the sky suggesting that the emission is 
dominated by the circumstellar envelope. We compare the ob- 
served azimuthally averaged intensity profiles of the MYSO 
envelopes to 1-D dust radiative transfer calculations. For each 
target, we attempt to reach consistency between the spatially re- 
solved 24.5 pm emission and the emission seen at other wave- 
lengths as represented by the SED. In addition, we describe the 
24.5 pm emission morphology in relation to known star forma- 
tion activity in the following subsection. 

3.2. Method of analysis 

Model SEDs and images are calculated using DUSTY, a code 
that solves the scaled 1-D dust radiative transfer problem 
(see Ivezic & Elitzur 1997). We use a spherically symmet- 
ric dust distribution illuminated by a central, unresolved star. 
Numerical solutions are independent of the star's bolometric 
luminosity, and the SED and images need to be scaled before 
making a comparison to the observations. The bolometric lu- 
minosity is the prime stellar parameter that sets the inner dust 
sublimation radius and thus the size scale of the envelope; an 
increase causes the size of the emitting region to be larger (as 
'"subi K V^boi)- As a result the intensity profile strongly de- 
pends on the Lboi assumed for the model. In practice, we de- 
termine the model Lboi by minimising the difference between 
the model and observed SED. The model Lboi can differ sig- 
nificantly from the observed L bo i, but for each model fit our 
adopted procedure finds the closest match between the two. In 
some case this match is poor and results in a large difference 
in observed and model Lboi- Observed L DO i are listed for each 
target in Col. 3 of Table 4. 

Other important model parameters are the outer radius of 
the envelope (R out ) and the total amount of dust. The latter is 



parametrised in DUSTY by the total optical depth Ay. Higher 
values for these two parameters tend to increase the long wave- 
length flux levels. Of course, neither the total amount of dust 
nor the envelope's outer radius can have any arbitrary value. 
We constrain the extinction by the 9.7/im silicate absorption 
feature that is generally found to be strongly in absorption 
among MYSOs. The required data were in most cases provided 
by mid-IR spectra taken with the short wavelength spectrom- 
eter (SWS, de Graauw et al. 1996) on board the ISO satellite 
(Kessler et al. 1996). The ISO data have been obtained from 
ESA's ISO data archive. 2 

The MYSO SEDs are constructed form the measured 
COMICS fluxes (Table 2) and from literature data. For most 
targets continuum measurements in the IR and (sub)mm are 
taken from the catalogue compiled by Gezari et al. (1999, 
available through Vizier), and includes IRAS and MSX ob- 
servations. These data were supplemented, where possible, 
with more recent observations, especially the compilation by 
Mueller et al. (2002). The continuum slope of the ISO-SWS 
spectrum longwards of the silicate absorption is used to com- 
plete the SED at mid-IR wavelengths. In three cases Spitzer 
MIPS data are available. We extracted the photometry apply- 
ing the non-linearity correction recipe by Dale et al. (2007). 
Data taken with large beams (> 15") were generally avoided in 
the model fitting procedure. The following subsections, that are 
devoted to each object, highlight the used and discarded data in 
the fitting process. In the accompanying figures for each source 
open symbols are used to indicate large beam data not taken 
into account in the model fitting procedure, whereas filled sym- 
bols indicate the data that were actually used in the fitting pro- 
cedure. Fluxes at wavelengths larger than 1.3 mm are excluded 
as they might be contaminated by free-free emission from ion- 
ized winds (e.g. Gibb & Hoare 2007). Overall, the constructed 
SEDs cover a wavelength range from about 1 pm to 1 mm. 

The model and observed images are compared in the fit- 
ting procedure by means of the azimuthally averaged intensity 
profile, normalized to its peak intensity. Model images are first 
scaled to model Lboi and convolved with the instrumental PSF. 
The intensity profiles are determined by binning and averag- 
ing the pixel counts in radial distance bins of size 0.13". The 
centre of the profile is found by minimising the pixel scatter 
as function of radial distance. We indicate the range in pixel 
counts found at each radius bin by an errorbar in the presented 
intensity profiles. 

In order to perform a systematic comparison between mod- 
els and observations, we have generated a standard grid that 
consists of 120000 DUSTY models. The grid probes the enve- 
lope parameter space and ranges in the following way: 

- Five radial density profiles. The radial density profile of the 
dust is described by a powerlaw of the form n(r) oc n r~ p . 
The five density profiles cover power indices from p = 0.0 
to p = 2.0 in steps of 0.5. 

- Four types of dust mixtures. The dust size distribution 
is MRN (Mathis, Rumple, & Nordsieck 1977), and re- 
mains untouched. We use DL (Draine & Lee 1984) opac- 

2 Accessible at this URL http://iso.esac.esa.int/ida/ 



6 



W.J. de Wit et al.: Subaru imaging of MYSOs 



ities for silicates and combine this with either DL opaci- 
ties for graphites or amorphous carbon (Zubko et al. 1996). 
Opacities are combined in a standard ISM mixture of 53% 
silicate and 47% graphite or a "twice ISM" mixture, i.e. 
67% silicate and 33% graphite. Different mixtures may 
accommodate the silicate absorption profile for given to- 
tal dust mass. The dust mixtures used have an opacity at 
1.3 mm of k between 0.3 - 0.4 g -1 cm 2 , similar to the val- 
ues of other types of dust in the literature, e.g. Ossenkopf 
& Henning (1994). 

- Two dust sublimation temperatures (r su bi) : 1000K and 
1500 K. The value for this parameter is the only unsealed 
parameter in DUSTY. The value of r su bi has a large ef- 
fect on the location of the dust sublimation radius but only 
a small effect on the normalized intensity profile. This is 
mainly due to a similar temperature dependence through 
the cloud with a very steep initial drop followed by a grad- 
ual decline. The similar temperature stratification leads to a 
similar normalized intensity profile. The different dust sub- 
limation temperatures produce a different balance between 
short and long wavelength flux in the SED. 

- Three scaled sizes for the envelope. The outer radius is 
scaled to the inner sublimation radius. The standard model 
grid contains models with sizes 750, 1000 and 1250 times 
the inner dust sublimation radius, which correspond to typ- 
ical envelope sizes of about 0.1 pc. 

For each combination of these four envelope structure pa- 
rameters, DUSTY calculates the resulting SEDs and images for 
an increasing total dust mass contained in the envelope. The Ay 
grid goes from 2 to 200 in steps of 2, resulting in a grand to- 
tal of 120 000 models. The envelope parameters adopted in this 
grid lead to typical temperatures found at the outer radius of 
20 K for r suM = 1 000 K for and 30 K for r subl = 1500 K with 
typical densities at R oM of 10 3 ~ 4 crrT 3 

The fitting procedure tries to find a model that fits the in- 
tensity profile and SED simultaneously. The SED fit includes 
the silicate absorption feature, slope of the ISO-SWS spectra, 
the far-IR SED peak and the (sub)mm continuum data. We 
aim to match the flux levels of the integrated (sub)mm con- 
tinuum data providing additional constraints on the model pa- 
rameters. The ISO spectra were rebinned logarithmically in 35 
wavelength bins between 2.5 and 43 yum in order to reduce its 
weight with respect to the continuum data. The binned spec- 
trum probes the silicate feature at 9.7 /mi and the change in 
continuum slope longwards of it. In short, the fit procedure ini- 
tially estimates the model Lboi by matching the overall shape 
of the scaled model SED to the observed one. DUSTY's output 
images are then accordingly scaled and convolved with the in- 
strumental PSF. A comparison of the model SED and intensity 
profile to the observed ones is made for all generated models. 
A simple tally is performed based on a goodness-of-fit crite- 
rion. We apply a least-squares criterion for the SED fit, and a 
weighted least-squares criterion for the intensity profile fit. The 
latter uses weights that are inversely proportional to the square 
of the range in normalized intensity at each radial distance bin. 
This range is represented by errorbars in the intensity profile 
plots for each MYSO target. The model that fits both sets of 



data best is the one with the highest average ranking in the two 
tallies. If no simultaneous satisfactory fits were obtained, we 
changed the emphasis of the procedure in order to get fits to in- 
tensity profile and the long wavelength range of the SED only, 
i.e. dropping the requirement of fitting the near-IR and mid-IR 
part of the SED. In some cases no satisfactory model was found 
within our standard grid, and we chose to explore the effects of 
changes in envelope size. This is discussed in more detail in the 
subsections dedicated to each source. 

Finally, the parameters for the underlying star are the same 
for each model and are the following. Lboi of the MYSOs are in- 
dicative of early B-type stars, thus we adopt a T e g of 25 000 K. 
There is only a negligible effect on the model dust emission for 
stars with slightly different effective temperatures. Distances 
for our sample stars were taken from the literature (see Table 4). 

All final model parameters are listed in the modelling 
overview Table 4. For each MYSO we give the parameters for 
three models. Models #1 correspond to the overall best-fitting 
model. Models #2 and #3 constitute the best-fit for models re- 
stricted to density powerlaws with indices that bracket the pow- 
erlaw index of the best-fitting model #1, and are not the second 
and third best-fitting model. This is to illustrate how well the 
powerlaw index is constrained by the observables as this is our 
main priority; in some cases, models #2 and #3 do not provide 
acceptable fits. Therefore, the model parameters of Table 4 can- 
not be considered to represent the uncertainties on the model 
parameters. More detailed explanation can be found in the sub- 
sections dedicated to each MYSO. 

3.3. 24.5 fim morphology and modelling 
3.3.1. S140(Figs.3and 4) 

Description: At a distance of 910pc, S140 is a bright-rimmed 
cloud that forms the interface between an H n region and the 
molecular cloud L1204 (Crampton & Fisher 1974). S140 con- 
tains an IR cluster of at least three sources (Beichmann et al. 
1979). The main CO outflow in the region goes SE-NW and a 
monopolar reflection nebula emanates from IRS 1 which is as- 
sociated with the blue-shifted SE lobe (Hayashi et al. 1987; see 
also Hoare & Franco 2007). The region displays spectacular 
arc-like features that are seen in high angular resolution TT-band 
images (Forrest & Shure 1986). They probably constitute ma- 
terial swept-up by outflow activity in the region (Weigelt et al. 
2002). Recent cm continuum and /f-band polarimetric imaging 
corroborate to suggest the presence of a disk-like structure for 
IRS 1 oriented in a NE-S W direction, perpendicular to the main 
outflow (Hoare 2006; Jiang et al. 2008). 
Mid-IR morphology: The COMICS image of the S140 region 
is presented in Fig. 3. It reveals a wealth of features associ- 
ated with previously known objects, the positions of which are 
labelled in the figure. Discrete peaks in mid-IR emission are 
found at the positions of IRS1, IRS2, and IRS3. The bright 
emission within 3-4" of IRS 1 has a fairly symmetric profile on 
the sky. The bright emission region 17" to the North of IRS1 
coincides closely, but not exactly with the position of IRS2N 
according to the positions given in Tofani et al. (1995). IRS2S 
is identified with a point-like mid-IR source just South of that. 
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Fig. 3. COMICS 24.5 yum image of the S140 region. The image is 
linearly scaled. Annotated objects are discussed in the text. Crosses 
correspond to radio sources (Evans et al. 1989; Tofani et al. 1995). 
Contour levels are at 0.5%, 1%, 2%, 3%, 5%, 10%, and 40% of peak 
flux density (5.9 10 2 Jy arcsec~ 2 ). North is up, East is to the left. 



IRS3 is found to be a triple system (Preibisch et al. 2003) and 
our image partly resolves the secondary object (IRS3b) at a dis- 
tance of 0.75" East of the primary source (IRS3a). 

The COMICS image also shows faint structures of diffuse 
emission. An arc of mid-IR emission concurs with the /T-band 
arc labelled "I" by Weigelt et al. (2002). The curved wisp at 
~ 3" from IRS3 corresponds to feature "F" of Preibisch et al. 
This feature is also known to have strongly polarised emis- 
sion and strong H2 line emission, suggesting scattering and 
shocked material. Patches of diffuse emission are also found 
coincident with the radio sources VLA4 and NW (see Evans et 
al. 1989). This is the first time that the two radio sources are 
seen in the mid-IR. Previously (Evans et al. 1989), they were 
found coincident with the brightest parts of very extended near- 
IR nebulosity. Finally, we note that the two submm emission 
peaks (Minchin et al. 1995; Thompson et al. 2006) do not have 
24.5 jt/m emission counterparts. Clearly the image reveals that 
24.5 jt/m emission is not restricted to compact MYSO envelope 
emission. The emission has a faint diffuse character when it is 
associated with shocked dense material. 

Model results: We focus the modelling on the main component 
in the region: bright MYSO IRS1. Its azimuthally averaged 
24.5 jt/m intensity profile is presented in panel a of Fig. 4. The 
errorbars cover the peak to peak range in pixel counts measured 
at each radial distance bin. Their small values indicate that 
IRS 1 can be considered to be symmetric to first order. We build 
the IRS1 SED (Fig. 4) from literature data. Photometry ob- 
tained using image restoration techniques applied to KAO con- 
tinuum observations at 50 and 100/im by Lester et al. (1986) 
are preferred over large beam data (50") presented by Schwartz 
et al. (1983) and IRAS. We adopt for IRS 1 the IRAM measure- 
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Fig. 4. Simultaneous model fits to the intensity profile and SED of 
S140IRS1. Panel a) Observed (with errorbars) and three model nor- 
malized azimuthal intensity profiles at 24.5 yum. The code for the 
model line styles is the same as in panel b. The errorbars indicate the 
range in observed intensity values found at each radial distance bin. 
The horizontal dashed line equals three times the rms level of the nor- 
malized background level. Panel b) Observed (symbols) and model 
SEDs (various line styles). The COMICS flux measurement is repre- 
sented by an asterisk. Open squares are IRAS measurements. Open 
triangles are large beam measurements discussed in the model results 
paragraph for each MYSO target. The thick line denotes the ISO-SWS 
spectrum. Models are fitted to the filled symbols and the logarithmi- 
cally binned ISO-SWS spectrum only. The best fit model (model #1 
in Table 4) is indicated by the full line, models #2, and #3 by a dashed 
line, dot-dashed linestyle respectively. The best fitting model (#1) has 
a p = 1.0 radial density profile. 



ment at 1.3 mm (11" HPBW) by Gurtler et al. (1991). Large 
beam far-IR and submm data by Schwartz et al. (1983) and 
Gurtler et al. (1991) are indicated by triangles in Fig. 4. 

The intensity profile and the SED can be reproduced simul- 
taneously by models with a p = 1.0 radial density distribu- 
tion. It is clear from Fig. 4 that the intensity profile can also 
be fit by model #2, and to a somewhat lesser extent, model 
#3. However, the simultaneous fit to the SED excludes shallow 
p = 0.5 (model #2) and steep p = 1.5 (model #3) models as 
viable alternatives. This is especially evident from their poor 
reproduction of the observed silicate absorption profile and the 
SED peak. On larger scales, modelling of the SED longwards 
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Fig. 5. COMICS 24.5 image of the M8E region. Relative posi- 
tions of the MYSO and the cometary shaped Hn region are consis- 
tent with near-IR (2MASS) and radio images (Simon et al. 1984). 
Contour levels are at 1%, 2%, 5%, 10%, and 40% of peak flux density 
(2.4 10 2 Jy arcsecr 2 ). North is up, East is to the left. 
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of 100 by Van der Tak et al. (2000) resulted in a p = 1.5 
density distribution. Mueller et al. (2002) find an overall best- 
fit to the SED and 350 fim intensity profile for spherical models Fig. 6. As Fig. 4, but for M8E-IR. The SED data is described in 



with p = 1.25. 



Sect. 3.3.2. The best fitting model has p = 1.25 radial density profile. 



3.3.2. M8E(Figs.5and 6) 

Description: M8E consists of a bright IR source and an opti- 
cally thin radio source 8" to its North- West (Simon et al. 1984). 
The region shows evidence of a bipolar CO outflow, but it is not 
clear which of the two sources powers it (Simon et al. 1984). 
The IR object is resolved on scales of milli arcseconds by lu- 
nar occultation observations at 3.8 fim and lOjjm. It can be in- 
terpreted as consisting of two physically different components 
(Simon et al. 1985): a compact hot component and a broad 
cooler component. The broad component, that has an angular 
size of 0.1", does not correspond to the 30 K molecular cloud 
material causing the SED to peak at 100 fim, but is suggested to 
be a transition region between relatively hot disk material and 
cold cloud material. 

Mid-IR morphology: Fig. 5 shows the detection of M8E-IR as 
a compact discrete source which is symmetric to first order. 
The radio source emission is diffuse and has a cometary mor- 
phology. Its peak emission is found at a distance of 7.7" from 
M8E-IR. A third source is detected 6" West of M8E-IR. 
Model results: The primary source of M8E-IR continuum data 
is from Mueller et al. (2002) and Gurtler et al. (1991). The six 
(sub)mm continuum measurements have been done with var- 
ious beamsizes between 19" and 30", except the 870//m one 
that was performed with a 9" half-power beamwidth (Gurtler 
et al. 1991). For reference, we also plot the 1.2 mm datapoint 



obtained by Beltran et al. (2006) with a half-power beamwidth 
of 26" (open square). It is clear that the intrinsic (sub)mm flux 
levels are uncertain, given the spread in flux levels in this wave- 
length range. At the short wavelength region we use the mea- 
surements by Simon et al. (1985). 

M8E-IR is compact and only marginally resolved. In fact, 
the intensity profile is too steep for the bolometric luminosity 
derived from the observed SED. Steeper radial density power- 
laws (creating a more compact cloud structure) lead to steeper 
intensity profiles, yet this advantage is offset by a decrease in 
flux at the long wavelength end. A commensurate increase in 
bolometric luminosity is required to fit the (sub)mm, but this in 
turn makes the intensity profile too shallow. Choosing a smaller 
outer radius is marred with the same problem, as is increasing 
the dust sublimation temperature. 

The model that fits the intensity profile best has a p = 1 .5 
powerlaw. Density profiles as shallow as p = 1 .0 are incompat- 
ible with the intensity profile. In order to increase the (sub)mm 
flux levels (without increasing the bolometric luminosity) we 
make the cloud larger than the sizes assumed in our model 
grid. Depending on the intrinsic (sub)mm flux levels, p = 1 .25 
to p = 1.5 are preferred with an outer radius of 3000 times 
the dust sublimation radius, or ~ 1 pc. The flux levels of the 
ISO-SWS spectrum cannot be reproduced by any model SED. 
On scales of 1000-10000 AUs, the 350//m intensity profile is 
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Fig. 7. COMICS 24.5 fim image of the AFGL2136 region. The large 
crosses correspond to the peak emission of the three near-IR scattering 
nebulae (Kastner et al. 1994). Contour levels are at 2%, 5%, 10%, and 
40% of peak flux density (1.1 10 2 Jy arcsec ). North is up, East is to 
the left. 

best represented by models with p = 1.75 radial density dis- 
tributions (Mueller et al. 2002), i.e. too steep to reproduce the 
24.5 pm intensity profile. 

3.3.3. AFGL 2136 (Figs. 7 and 8) 

Description: A three-lobed near-IR reflection nebula littered 
with some 30 faint point sources is illuminated by the 
dominant source IRS1. Near-IR polarisation indicates that 
IRS1 is illuminating conical cavities within a dusty enve- 
lope (Kastner, Weintraub, & Aspin 1992). IRS1 is the driving 
source of an arcminute-scale bipolar CO outflow with a R A. of 
135° (Kastner et al. 1994). Weak, optically thick radio emission 
originating from IRS1 is detected by Menten & van der Tak 
(2004). The radio source is somewhat elongated in the South- 
Easterly direction. 

Mid-IR morphology: The COMICS image in Fig. 7 reveals one 
dominant source, IRS1, that is resolved and symmetric. The 
image is not very deep, and the envelope can be traced out to 
~3", where noise starts to dominate. IRS1 displays an extended 
wing to the South-East coincident with the near-IR "South" 
lobe of the nebula. A faint counter wing to the North- West 
would correspond to the "West" lobe. There are also traces of a 
24.5 pm counterpart to the faint "East" lobe (terminology from 
Kastner et al. 1992). The bright South lobe extends at roughly 
the same P.A. as the outflow activity. 

Model results: The SED is built using the set of (sub)mm data 
from Kastner et al. (1994) obtained with the JCMT (16.8"- 
18.5" HPBW), except the datapoint at 350 pm which is from 
Mueller et al. (2002). IRAS measurements agree closely with 
the Harvey et al. (2000) data at 50 and lOOyum . 
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Fig. 8. As Fig. 4, but for AFGL 2136. The SED data is described in 
Sect. 3.3.3. The best fitting model has p = 1.0 radial density profile. 

The intensity profile and SED of AFGL 2 136 are best fit 
by a p = 1.0 model. None of the models that can reproduce 
the intensity profile do a particular good job in reproducing 
the SED's short wavelength range, i.e. < 10 pm. At larger spa- 
tial scales, Mueller et al. (2002) and van der Tak et al. (2000) 
prefer steeper models with p = 1.75 and p = 1.25, respec- 
tively. Harvey et al. (2000) find reasonable fits to the 50 and 
100/im intensity profiles and the SED, by dividing a p = 1.5 
model envelope in a high and low optical depth "hemisphere". 
We also find that certain p = 1 .5 models are capable of nicely 
fitting the full SED longwards of 5 pm (not shown in Fig. 8), but 
they do not reproduce the 24.5 pm intensity profile. The p = 1 .5 
model shown in Fig. 8 does reproduce the intensity profile, but 
fits the SED worse than our preferred p = 1.0 model. 

3.3.4. AFGL2591 (Figs. 9 and 10) 

Description: Various aspects of the AFGL 2591 envelope are 
discussed in van der Tak et al. (1999) using an extensive set 
of line and continuum observations. The authors conclude that 
the mid-IR and far-IR emission originates from the envelope, 
and that the line-of-sight to the source nearly coincides with 
the opening cone carved by an East- West oriented bipolar CO 
outflow. There is ample evidence that the circumstellar material 
is not distributed in a spherically symmetric fashion and that a 
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Fig. 9. COMICS 24.5 pm image of the AFGL2591 region. The rela- 
tive positions of the MYSO and the two cometary shaped H n regions 
are consistent with the 8 GHz sources in Tofani et al. (1995). Contour 
levels are at 0.5%, 1%, 2%, 5%, 10%, and 40% of peak flux density 
(7.7 10 2 Jy arcsec~ 2 ). North is up, East is to the left. 

low-opacity pathway close to the line of sight has important ef- 
fects on the source's appearance (Preibisch et al. 2003). 
Mid-IR morphology: The COMICS image presents three 
sources: bright IRS1, a cometary shaped emission feature to 
the South- West, and a small emission region to the North- West. 
IRS 1 is so bright that the spider diffraction pattern and the de- 
tector cross-talk (the vertical dark lane) are visible in Fig. 9. 
The cometary shaped feature to its South-West closely fol- 
lows the optically thin H n region VLA1 (Wynn- Williams et al. 
1977; Trinidad et al. 2003). The relative position of the source 
at the Northern rim with respect to the other two is consistent 
with 8 GHz source "n4" from Tofani et al. (1995). 
Model results: In Fig. 10, we use JCMT submm data from 
Jenness et al. (1995), and CSO data from Mueller et al. (2002), 
that have comparable beamsizes (18" and 14", respectively). 
At mm wavelengths we show data from Giirtler et al. (1991) 
and one datapoint at 1.3 mm from Walker, Adams & Lada 
(1990). The Giirtler et al. data were taken with a beam of 
~10" whereas the Walker et al. data have a 30" beam, possibly 
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Fig. 10. As Fig. 4, but for AFGL 2591 IRS 1. The SED data is de- 
scribed in Sect. 3.3.4. No simultaneous model fit to SED and intensity 
profile is found. 



explaining the flux level difference between these two sets. The 
far-IR KAO data of Lada et al. (1984) closely corresponds to 
the IRAS 60 fim and 100 /mi photometry. 

AFGL 2591 IRS1 presents an interesting case, as no si- 
multaneous SED and intensity profile fit can be obtained. The 
bolometric luminosity implied by the observed SED impedes 
any reasonable model fit to the intensity profile (see e.g. the 
dot-dashed p = 1.0 model in Fig. 10a and b). We attempt to 
find a solution by simply fitting the (sub)mm SED and inten- 
sity profile. The compromise that does fit the 24.5 pm flux and 
(nearly) the intensity profile is given by a p = 1 .5 model with 
a large outer radius of ~ 5' (1.6pc). The discrepancy found be- 
tween spherical models and AFGL 2591 is not limited to the 
mid and near-IR but is still severe at far-IR wavelengths. Such 
models require a substantial revision of the object's bolomet- 
ric luminosity. In other words, the source is more compact than 
predicted by spherical models. A demonstration of a more ap- 
propriate 2-D modelling applied to this source is reported in 
Preibisch et al. (2003). These authors using speckle interfer- 
ometry probing scales of 170 AU. They are able to model a 
resolved structure either as the inner rim of a circumstellar disk 
or as the dust sublimation radius of the MYSO envelope. 
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Fig. 11. COMICS 24.5 yum image of NGC2264 IRS1. Contour lev- 
els are at 1%, 2%, 5%, 10%, and 40% of peak flux density 
(4. 1 10 2 Jy arcsecT 2 ). North is up, East is to the left. 



3.3.5. NGC2264IRS1 (AFGL989, Figs. 11 and 12) 

Description: A multi-wavelength study by Schreyer et al. 
(2003) suggests that NGC 2264 IRS 1 is a young B-type star 
surrounded by low-mass companions located in a low-density 
cavity of a clumpy, shell-like, and dense cloud remnant. 
NGC 2264 IRS 1 shows evidence for a CO outflow oriented 
along our line of sight (Schreyer et al. 1997). Attempts to re- 
solve sub-arcsecond structure related to the CO outflow cav- 
ities in the near-IR with HST (Thompson et al. 1998) and us- 
ing speckle techniques (Alvarez et al. 2004) failed. A reflection 
nebula on scales of arcseconds is clearly visible in the optical 
(Scarrott & Warren-Smith 1989) and near-IR (Schreyer et al. 
1997). 

Mid-IR morphology: The COMICS image shows a single, sym- 
metric object (IRS1). Most of the substructure seen in the im- 
age are diffraction patterns. 

Model results: The model fits to SED and intensity profile 
are presented in Fig. 12. IRS1 coincides with submm source 
MMS5 (Ward-Thompson et al. 2000; Peretto et al. 2006) 
if one takes the correct near-IR coordinates (from 2MASS: 
a = 06 h 41 m 10.16 s and 5 = +09°29'33.7") and the phase- 
referenced map by Nakano et al. (2003). We emphasise that 
the often reported lack of (sub)mm emission of IRS1 is due 
to wrong coordinates for IRS1 combined with poor astrome- 
try (Ward-Thompson et al. 2000; Schreyer et al. 2003; Nakano 
et al. 2003; Peretto et al. 2006). Far-IR (25 /mi and onwards) 
observations presented in the literature are done with beams 
larger than 20" (Harvey, Campbell & Hoffmann 1977; Chini 
et al. 1986; IRAS) and are considered to be upper limits. The 
model fit uses the JCMT and IRAM (sub)mm data presented 
in Ward-Thompson et al. (2000) with beam FWHM of 6", 8", 
13" and 12" for the 350, 450, 800 and 1300 pm data, respec- 
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Fig. 12. As Fig. 4, but for NGC 2264 IRS1. The SED data is described 
in Sect. 3.3.5. The best fitting model has p = 1.5 radial density profile. 



tively. Archive Spitzer 70 pm MIPS data are available for IRS 1 . 
We extracted the photometry applying the non-linearity correc- 
tion recipe of Dale et al. (2007). A 70yum flux value of 960 Jy 
for IRS1 is used in the model fits. 

The modelling indicates that bolometric luminosities that 
generate reasonable model fits to the intensity profile, find a 
rough correspondence with the mid-IR flux levels only if the 
powerlaw exponent equals p = 1.5. The 10pm MIPS is not 
attained by any model. It is clear that any spherical model 
aimed at fitting the far-IR or (sub)mm part of the SED would 
need much larger luminosities, which is incompatible with the 
24.5 pm intensity profile. 

3.3.6. S255 (Figs. 13 and 14) 

Description: The region is dominated by two near-IR sources 
(NIRS1 and 3) and their bipolar IR reflection nebulae that 
were initially identified as the single source IRS1 (Beichman 
et al. 1979; Tamura et al. 1991; Itoh et al. 2001). High resolu- 
tion near-IR speckle observations by Alvarez & Hoare (2004) 
resolve the bipolar reflection nebula of NIRS1 (the Western 
source), revealing it to be twisted in an "S" shape. High resolu- 
tion Subaru near-IR polarimetric data indicate a depolarisation 
plane perpendicular to the near-IR bipolar nebula, suggesting 
the presence of a disk (Jiang et al. 2008). 
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Fig. 13. COMICS 24.5 fim image of the S255 region. Contour lev- 
els are at 1%, 2%, 5%, 10%, and 40% of peak flux density 
(1.9 10 2 Jy arcsecT 2 ). North is up, East is to the left. 



Mid-IR morphology: The COMICS image shows NIRS3 to be 
the dominant source at 24.5 pm. Both mid-IR sources are found 
to be symmetric and resolved. 

Model results: We focus our analysis on the bright source 
NIRS3 (the Eastern source) for the SED and intensity profile 
fit. The spectral range of the ISO spectrum does not include 
the 9.7 pm silicate feature. Instead we use the spectrum pub- 
lished in Willner et al. (1982), that describes the depth of the 
feature, but does not constrain the slope of the continuum long- 
wards of it. The short wavelength range of the Willner et al. 
spectrum corresponds well with photometry taken by Evans et 
al. (1977). The SED uses (sub)mm data from the combined 
"core-envelope" emission at 350 pm with a FWHP of 30"from 
Metzger et al. (1988), and measurements (FWHP of 40") by 
Richardson et al. (1985). The mm measurement is by Chini et 
al. (1986a) obtained with the IRTF with a 90" beamsize. Other 
sets of (sub)mm points with higher fluxes are available in the 
literature and are represented by open triangles (Richardson 
et al. 1985; Metzger et al. 1988; Klein et al. 2005). The total 
COMICS flux level (i.e. NIRS1 and 3 taken together) is within 
5% of the MSX flux. More emission structure is hidden in the 
IRAS beam (see e.g. Longmore et al. (2006) for larger scale 
18.7 pm images). 

S255 IRS3 constitutes a source that is mildly resolved at 
24.5yum . The intensity profile for IRS3 has been obtained 
by excluding the region between P. A. -45°and -135°, where 
emission from IRS1 dominates. In Fig. 14 the observations are 
compared to three spherical envelope models, the parameters of 
which are listed in Table 4. The model with a radial density dis- 
tribution with power exponent p = 1.25 best fits both the SED 
and intensity profile. The p = 1.5 model is clearly excluded 
as it requires high optical depths not substantiated by the sili- 
cate absorption profile. The case of S255 IRS3 is similar to the 
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Fig. 14. As Fig. 4, but for S255 IRS3. The SED data is described in 
Sect. 3.3.6. The best fitting model has p = 1.25 radial density profile. 



one encountered for M8E-IR. The models are required to have 
large envelopes which increase the long wavelength flux. This 
can be achieved without changing the bolometric flux, keeping 
the intensity profile unchanged. 

3.3.7. AFGL5180 (Figs. 15 and 16) 

Description: AFGL 5180 has a blueshifted CO flow with a PA. 
of ~ 130° (Snell et al. 1988). Saito et al. (2006) show that the 
IRAS point source breaks up into two main (sub)mm cores (see 
also Minier et al. 2005). One of these cores is centred on the 
near-IR source NIRS1 as identified by Tamura et al. (1991). 
The region reveals an intricate collection of mm cores mixed in 
with various near-IR sources (Tamura et al. 1991; Saito et al. 
2007). 

Mid-IR morphology: At 24.5 pm and at the resolution of 
Subaru, AFGL 5 180 consists of three sources: one extended 
source and two point-like sources. We follow the nomenclature 
by Tamura et al. (1991), in which the main discrete source can 
be identified as NIRS1 and the diffuse source as NIRS2. The 
region was recently imaged by Longmore et al. (2006) in the 
mid-IR. Their deconvolved images at 7.9 pm shows that NIRS 1 
itself consists of three components that have a mutual distance 
of about 1". The COMICS image shows that NIRS1 is some- 
what asymmetric, but it is not resolved in multiple sources. 
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Fig. 15. COMICS 24.5 jum image of the AFGL 5 1 80 region. The prin- 
cipal Western source corresponds to NIRS1 (Tamura et al. 1991). 
Contour levels are at 5%, 10%, and 40% of peak flux density 
(4.5 10 2 Jy arcsec~ 2 ). North is up, East is to the left. 



Model results: We concentrate our discussion on NIRS1, the 
main mid-IR source. Model comparisons to the observations 
are less well constrained than in the previous cases, because 
NIRS1 lacks a 10 pm spectrum. An independent measure of 
the total optical depth is therefore absent. We use the (sub)mm 
photometry from Gear et al. (1988) and Thompson et al. 
(2006). The Gear et al. data have FWHM beamsizes of about 
an arcminute, whereas the Thompson et al. data are taken 
with SCUBA with FWHM beamsizes of 8" and 14" for the 
450yum and 850yimi observations. The far-IR datapoint is from 
Ghosh et al. (2000). For the mid-IR we use the MSX photome- 
try, that have a good correspondence with the IRAS data points. 
We find reasonable fits to the intensity profile for p = 0.0, 
p = 0.5 and p = 1.0 radial density distributions (see Fig. 16), 
but the MSX and IRAS points in the SED are better reproduced 
by the p = 1 .0 model. Any steeper distributions may provide 
equally good fits, except that they would require optical depths 
exceeding Ay = 200 m and much larger outer radii then adopted 
in our grid, which are typical for most sources in our sample. 
The p = 1 .0 model already requires a comparatively large op- 
tical depth with respect to the other MYSO envelopes for it to 
fit the MSX data. 

3.3.8. IRAS 201 26+41 04 (Figs. 17 and 18) 

Description: IRAS 20126 is one of the best examples of a rel- 
atively massive young star that is actively fed by an accretion 
disk. The star is thought to have a mass of approximately 7 M 
(Cesaroni et al. 2005). The system is known to show jet/outflow 
phenomena at a RA. of about -60° (see Cesaroni et al. 2005; 
Su et al. 2007), with the jet oriented almost perpendicular to the 
line of sight (Cesaroni et al. 1999). At wavelengths shorter than 
20 pm a "dark lane" separates two emission regions (Sridharan 
et al. 2005; De Buizer 2007). This morphology could be either 
due to the presence of two tight clusters (in which case the dark 
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Fig. 16. As Fig. 4, but for AFGL 5 180. The SED data is described in 
Sect. 3.3.7. The best fitting model has p = 1.0 radial density profile. 



lane does not correspond to a physical structure) or due to a sil- 
houetted disk with emission emanating from the cavity walls 
of the outflow either side of the disk (De Buizer 2007). 
Mid-IR morphology: The COMICS image consists of two 
emission regions. The morphology of the North- West source 
(source "5" in De Buizer 2007) corresponds to the one pre- 
sented in Shepherd et al. (2000) at 17.9/Lan and coincides with 
the direction of the outflow. The South-Eastern peak emission 
overlaps the two sources that are seen separated by a dark lane 
at shorter wavelengths. No separation due to a dark lane is de- 
tected, although it is extended perpendicular to the dark lane 
seen at shorter wavelengths. The two sources have a separation 
of 0.7" in 18.3 pm images by De Buizer (2007). The somewhat 
lower resolution at 24.5 micron is probably partly responsible 
for the failure to separate the two sources. In addition, it seems 
likely that lower extinction and thermal emission from the pu- 
tative dark lane itself could be blurring the distinction between 
the two emission components. It would argue in favour of the 
dark lane being a physical structure rather than a clearance of 
emitting material as suggested in De Buizer (2007). 
Model results: The ISO-S WS spectrum has a low SNR at wave- 
lengths shorter than 10 pm. This includes the silicate absorption 
profile, inhibiting a independent measure of the envelope's op- 
tical depth. A sharp rise of the spectrum towards longer wave- 
lengths indicates a very high extinction. In the model fits we 



14 



W.J. de Wit et al.: Subaru imaging of MYSOs 



... 



■ ■ ■ ■ 




1 -1 

Arc Seconds 



Fig. 17. COMICS 24.5 pm image of IRAS 20126+4104. Contour lev- 
els are at 5%, 10%, 20%, 40%, and 70% of peak flux density 
(4.7 10 1 Jy arcsecT 2 ). North is up, East is to the left. 



use the photometry collected by Hofner et al. (2007, and ref- 
erences therein). We fit only the large beam data presumably 
corresponding to the dusty halo. The submm observations are 
obtained with the JCMT by Cesaroni et al. (1999) with HPBW 
between 7" and 14". The intensity profile is obtained from the 
region that excludes emission from the North- Western source. 
We find good simultaneous fits to both intensity profile and 
SED, with a preference for flat radial density distributions with 
powerlaw indices p = 0.0 and p = 0.5. Steeper density distri- 
butions are not able to fit the 24.5 pm intensity profile. These 
values for the density distribution power index are at odds with 
Hofner et al. (2007) and van der Tak et al. (2000). The for- 
mer find consistency between the observed SED and what one 
would expect from an accretion disk embedded in a spherical 
infalling halo (on scales of 10") with p = 1.5. Van der Tak et 
al. (2000) find p = 1.75, on scales larger than the ones probed 
by the 24.5 pm image. 



3.3.9. MonR2(Figs.19and 20) 

Description: The central region of the Mon R2 cloud consists 
of at least five bright IR sources within ~ 0.25 pc (Beckwith 
et al. 1976). The region shows complex submm dust contin- 
uum emission which is well illustrated by the 850 pm map by 
Giannakopoulou et al. (1997). High resolution near-IR imag- 
ing (0.075" resolution) resolves the brightest mid-IR source 
(IRS3) into a triple system surrounded by strong diffuse neb- 
ulosity (Preibisch et al. 2002). 

Mid-IR morphology: The COMICS mosaic in Fig. 19 shows 
IRS3 to be the dominant discrete source at 24.5 pm. The large, 
30" scale mid-IR shell structure is identifiable in near-IR con- 
tinuum emission in which case it is due to dust scattering 
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Fig. 18. As Fig. 4, but for IRAS 20126+4104. The SED data is de- 
scribed in Sect. 3.3.8. The best fitting models have p = 0.0 or p = 0.5 
radial density profiles. 



(Howard et al. 1994). Its morphology corresponds roughly to 
the extent of the blister Hn region (Massi et al. 1985). If we 
compare the mid-IR image with these maps, we may iden- 
tify the 24.5 pm emission with the walls of the ionized region. 
Maximum intensity along the mid-IR ridge corresponds to the 
location of the source IRS1. The image does not reveal any 
point source at this position. IRS2, IRS3 and IRS5 are all ob- 
served to be extended sources. IRS3 is elongated with a P.A. 
along the direction of the main binary (Preibisch et al. 2002). 
Model results: Submm and mm imaging at 870jt/m and 1.3 mm 
(HPBW of 18" and 23", respectively) by Henning et al. (1992) 
resolve and identify the various components of the cloud with 
the IR sources. We show in Fig. 20 the intensity profile and 
SED fits to IRS3, the most luminous mid-IR source of the re- 
gion. Fig. 20b also shows the Mueller et al. (2002) compila- 
tion of continuum measurements at the long wavelengths, in- 
cluding the 350/mi CSO measurement (HPBW 14"). The flux 
level of the ISO spectrum corresponds closely to the COMICS 
flux measurement. The intensity profile levels out at distances 
> 5" which is probably due to contribution in flux by the dif- 
fuse emission in the region. 

The case of Mon R2 IRS3 is similar to S140 IRS1. 
Simultaneous model fit to the intensity profile and SED shows 
that models with a radial density profile of p = 1.0 are pre- 
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Fig. 19. COMICS 24.5 yum image of the Mon R2 region. Contour 
levels are at 1%, 3%, 20%, and 70% of peak flux density 
(5.0 10 2 Jy arcsecT 2 ). North is up, East is to the left. 



ferred. Again steep radial profiles are excluded as they require 
high optical depths in order to provide a good fit to the intensity 
profile. 

3.3.1 0. AFGL 437S (Figs. 21 and 22) 

Description: AFGL 437 is a compact cluster of few dozen near- 
IR sources located at a distance of 2.7 kpc (Wynn- Williams 
et al. 1981; Weintraub & Kastner 1996). The cluster is dom- 
inated by four bright sources 437N, S, E and W. The source 
437W dominates the radio emission, whereas only weak radio 
emission is measured towards the 437S. (Kurtz, Churchwell, & 
Wood 1994; Torrelles et al. 1992). Alvarez et al. (2004) resolve 
a monopolar sub-arcsecond near-IR nebula from 437S. Water 
masers have been detected towards 437W and 437N (Torrelles 
et al. 1992). Weintraub & Kastner (1996) find that 437N ac- 
tually breaks up into two components (see also Rayner & 
McClean 1987). The South-Eastern source of the two (named 
WK34) is the most embedded and found to be responsible for 
the region's molecular outflow and the dominant source of the 
near-IR reflection nebula (see also Meakin et al. 2005). 
Mid-IR morphology: The COMICS image is dominated by the 
UCH n region associated with 437 W. No discrete counterpart 
is found for 437W, only a diffuse emission region. The 437N 
source is resolved into two components, with the driving source 
WK34 the most luminous of the two. Of the 4 sources detected, 
437S is the most luminous and marginally resolved by our ob- 
servations at 24.5 pm. 

Model results: We concentrate our analysis on 437S. Submm 
and mm observations were performed by Dent et al. (1998) 
with the JCMT at four wavelengths with beamsizes between 
16" and 19". The ISO-SWS spectrum of the region is domi- 
nated by emission from the H n region and is discarded. Instead 
we use the 10pm photometry of 437S taken by Wynn- Williams 
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Fig. 20. As Fig. 4, but for Mon R2 IRS3. The SED data is described in 
Sect. 3.3.9. The best fitting model has p = 1.0 radial density profile. 



et al. (1981). The model fitting procedure discards the IRAS 
photometry, but takes into account the COMICS flux measure- 
ment. 

The spatial information derived from the marginally re- 
solved source is not enough to strongly constrain the various 
radial density distributions. We show in Fig. 22 three models 
that reproduce the intensity profile and the continuum emis- 
sion in the SED. We find that p = 1.0 profiles require large 
optical depth in order to fit the (sub)mm and the 24.5 pm data. 
This seems to be excluded by the silicate absorption profile. We 
therefore prefer rather shallow radial density distributions with 
p = 0.0 or p — 0.5, but with considerable uncertainty. 

3.4. Complex sources 

MYSO sources presented in this section show evidence 
for multiple condensations within 1" of the profile centre 
(AFGL 4029, AFGL 961, and W3IRS5). Azimuthally aver- 
aged intensity profiles yield consequently large uncertainties, 
which constrain the models only little. In the case of the fourth 
source presented in this section, the Cep A star forming region, 
no discrete central source could be identified at all. We discuss 
the 24.5 pm images alongside some literature background for 
each of these four objects. 
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Fig. 21. COMICS 24.5 yum image of the AFGL437 region. Contour 
levels are at 5%, 10%, 20%, and 50% of peak flux density 
(3.2 10 1 Jy arcsecT 2 ). North is up, East is to the left. 



3.4.1 . AFGL4029 (Figs. 23 and 24) 

Description: Beichman et al. (1979) have shown that the re- 
gion consists of two sources at 10 and 20 pm, separated by 
about 20". Both mid-IR sources have radio emission (Kurtz 
et al. 1994). IRS1 shows evidence for outflow activity. Zapata 
et al. (2001) find IRS1 to consist of a double radio source, in- 
terpreted to be a binary object. They identify the Southern of 
the two sources responsible for the outflow activity. 
Mid-IR morphology: The COMICS image reveals two emis- 
sion regions. A compact emission region identified as IRS1 
which is shown in Fig. 23, and a diffuse emission region as- 
sociated with IRS2 located at 20" from IRS1 (not shown 
in Fig. 23). IRS1 shows a discrete source with evidence for 
at least two condensations, the Southern being the bright- 
est one. In addition at least two bands of emission are lo- 
cated within the inner 2", rendering IRS 1 a rather patchy 
source at 24.5 pm. Finally, IRS 1 has an extended wing (~4") 
to the West. Scattered light emission and thermal emission 
due to extended nebulosity to the West of IRS1 is discussed 
in Deharveng et al. (1997) and Zavagno et al. (1999). The 
24.5 pm morphology of IRS1 follows closely the radio mor- 
phology as presented in Zapata et al. (2001). The mid-IR wing 
is possibly the counterpart to the East- West extension seen in 
high-resolution radio maps. 

Model results: The mid-IR part of the SED is adopted from 
data presented in Zavagno et al. (1999). IRAM 1.27 mm flux 
(beamsize ~15") is taken from Klein et al. (2005), and submm 
SCUBA fluxes (beamsizes of 8" and 14" for the 450yum and 
850yum) are presented in Di Francesco et al. (2008) . Finally a 
Spitzer MIPS flux measurement gives a flux density of 680 Jy 
for this source. 
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Fig. 22. As Fig. 4, but for AFGL437S. The SED data is described in 
Sect. 3.3.10. The best fitting model has p = 0.0 radial density profile. 



AFGL4029 presents an intensity profile that is character- 
ized by a large scatter in intensity, as represented by the error- 
bars in Fig. 24, due to the presence of patchy emission within 
2". The lack of a 10pm spectrum for the source also compli- 
cates the analysis. The data of Zavagno et al. (1999) cover the 
wings of the silicate absorption profile, providing some con- 
straints. Although the intensity profile shows a large scatter, 
simultaneous model fits to SED and intensity profile are hard 
to find; the basic problem being the extent of the source. This 
translates into a very luminous object, incompatible with the 
SED. In Fig. 24, we show models with relatively flat density 
profiles (p = 0.0 or 0.5), steeper profiles are incompatible with 
the intensity profile. 

3.4.2. AFGL961 (Figs. 25 and 26) 

Description: The source is a well-known double object located 
in the outskirts of the Rosette Nebula (Lenzen et al. 1984) and 
surrounded by a stellar cluster (Aspin 1998). The East and West 
components have a separation of about 6"and display outflow 
phenomena (Lada & Gautier 1982). Snell and Bally (1986) de- 
termined the radio spectral index, finding it to be consistent 
with an ionized wind. 

Mid-IR morphology: The dominant source in the COMICS im- 
age is AFGL961E. The Western source 96 1W is discrete and 
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Fig. 23. COMICS 24.5 jum image of the AFGL4029 region. Contour 
levels are at 5%, 10%, 25%, 50%, and 80% of peak flux density 
(2.2 10 1 Jy arcsecT 2 ). IRS2 is 20" to the East to IRS1, and not shown 
in the image. North is up, East is to the left. 

somewhat extended in a North- Westerly direction. AFGL 96 1 E 
displays two diffuse blobs to the North- West and South- West of 
the peak emissions. 

Model results: Chini et al. (1986) model this source and find 
good correspondence with a spherical model for wavelengths 
shorter than 100/im. An emission excess with respect to their 
spherical model is found for wavelengths longer than this. We 
use SCUBA measurements at 450 //m and 850 //m (beamsizes 
of 8" and 14"), which are significantly lower. Further, we 
use the spectrum taken by Willner et al. (1982), mid-IR pho- 
tometry from Cohen (1973), Simon & Dyck (1977), and the 
350/rni IRTF measurement (HPBW 90") is from Gurtler et 
al. (1991). We measured the 70/im flux from archival Spitzer 
MIPS data, finding a flux of 720 Jy for the source. The broad 
SED and intensity profile are well fit by p = 0.5 density pow- 
erlaws with a moderate optical depth. Steeper density profiles 
do not fit the intensity profile, and would additionally require 
either a higher optical depth or a much larger envelope outer 
radius in order to fit the submm points. 

3.4.3. W3(Fig.27 and 28) 

Description: W3 is an important and complex region contain- 
ing objects in various stages of the formation process (Wynn- 
Williams et al. 1972; Claussen et al. 1994). One prominent 
source in the region is IRS5, which consists of a double source 
(1" separation) at infrared wavelengths (Howell et al. 1981; 
van der Tak et al. 2005). IRS5 was discovered to harbour at 
least seven very compact radio continuum sources within a ra- 
dius of 0.03 pc (Claussen et al. 1994). A number of these are 
likely to be shock excited clumps in the surrounding molecular 
material, driven by embedded OB stars (Wilson et al. 2003 ;van 
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Fig. 24. As Fig. 4, but for AFGL 4029. The SED data is described in 
Sect. 3.4.1. The best fitting model has a p = 0.0 radial density profile. 



der Tak et al. 2005). Water masers in IRS5 have been shown 
to trace two outflows, both in roughly a North-South direction 
(Imai et al. 2000), a direction similar to the overall CO outflow 
of the region. 

Mid-IR morphology: The main 24.5 pm component in W3 re- 
gion is the bright source IRS5. The double source at its cen- 
tre is resolved, the Northern one is the brightest of the two 
sources. IRS5 shows prominent diffuse emission that surrounds 
the double source. It extends in both North-South and East- 
West direction. IRS6 is a resolved source with the peak emis- 
sion somewhat offset to the West of centre. IRS7 is a diffuse, 
cometary shaped emission region. IRS3 is resolved in various 
24.5 yum patches of diffuse emission. An unidentified source is 
found 10" North of IRS5. 

Model results: Data for W3 IRS5 are from van der Tak et al. 
(2005), the compilation made by Campbell et al. (1995), and 
the 10pm spectrum presented in Willner et al. (1982). The 
(sub)mm data have been obtained with beamsizes between 
15" and 19". 

Shallow radial density models provide good fits to the SED, 
but fail to reproduce the intensity profile. Models with p = 1.5 
powerlaw produce to much mid-IR flux requiring a high optical 
depth, leading to a too deep silicate absorption feature. The best 
fitting models require a p = 1.0 powerlaw. 
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Fig. 25. COMICS 24.5 /im image of the AFGL961 region. Contour 
levels are at 1%, 4%, 10%, 20%, 25%, and 50% of peak flux density 
(1.3 10 2 Jy arcsecT 2 ). North is up, East is to the left. 

3.4.4. CepA(Fig.29) 

Description: Patel et al. (2005) reported a flattened HCN struc- 
ture to be a circumstellar disk near the HW2 radio jet, the 
brightest radio source of the region (Hughes & Wouterloot 
1984). The radio jet is known to be the driving force of a large- 
scale molecular outflow in a North-East, South- West direction 
(Gomez et al. 1999). Jimenez-Serra et al. (2007) showed that 
the HCN disk is resolved into a much smaller disk and a hot 
core (Martin-Pintado et al. 2005). The centre of their disk falls 
close to the centre of the ionized HW2 jet as determined in 
Curiel et al. (2006). 

Mid-IR morphology: The region shows a complex morphology 
comprising various arcs and patches. The bright 24.5//m emis- 
sion corresponds to the bright near-IR reflection nebula that 
occupies the blue-shifted outflow cavity (Lenzen et al. 1984). 
The COMICS image does not show any point source and we 
performed a cross matching with a lower angular resolution, 
archival Spitzer 8 pm image in order to determine a rough as- 
trometric solution. The result shown in Fig. 29 is a reasonable 
match between 8 and 24.5 pm. Note that the 8 pm emission cor- 
responding to the 24.5 /jm emission peaks is saturated. HW2 
appears to be hidden by high extinction just at the bottom tip 
of a dark lane. 

Model results: Cep A is not modelled because direct radiation 
from the source is obscured at 24.5 pm. 

4. Summary of the model fitting 

We give a summary of the best-fitting models in Table 5. 
Whether the best model reproduces certain important SED 
wavelength intervals and the 24.5 pm intensity profile has been 
indicated by tickmarks and crosses. In some cases no such as- 
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Fig. 26. As Fig. 4, but for AFGL961E. The SED data is described in 
Sect. 3.4.2. The best fitting model has a p = 0.5 radial density profile. 



sessment could be made due to lack of data or data of insuffi- 
cient quality. The table makes clear that spherical models are 
capable of reproducing at least part of the SED and intensity 
profile simultaneously. The models systematically fail to repro- 
duce the short wavelength range, which is a well-known short- 
coming of spherical geometries (e.g. Giirtler et al. 1991). 

A pattern emerges from Table 5 in which the objects that 
are satisfactorily reproduced by spherical models (S140 IRS1, 
AFGL2136, AFGL5180, MonR2IRS3) are described by ra- 
dial density distribution with a p = 1 .0 power index. Four ob- 
jects that require steeper (p > 1.0) density distributions, viz. 
AFGL2591, NGC 2264 IRS 1, M8E-IR, S255IRS3, fail in re- 
producing either the silicate absorption or the intensity profile. 
These objects seem to show excess mid-IR flux with respect 
to spherical models that can only be accounted for if higher 
bolometric luminosities are adopted, a solution denied by the 
24.5 pm intensity profile. The MYSO that poses the biggest 
problem is AFGL2591. Van der Tak et al. (2000) forwarded 
the idea that the star's outflow activity is directed nearly along 
the line-of-sight, and the problems with spherical models en- 
countered here could be a manifestation of this. Schreyer et 
al. (1997) have argued for NGC 2264 IRS1 that again there 
is outflow activity along the line-of-sight towards this source. 
M8E-IR shows evidence of a CO outflow with the blue and red 
shifted components projected on top of each other, arguing for 
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Fig. 27. COMICS 24.5 pm mosaic of the W3 region. Contour lev- 
els are at 1%, 2%, 5%, 15%, 50%, and 85% of peak flux density 
(3.8 10 2 Jy arcsecT 2 ). North is up, East is to the left. 

a small inclination angle (Simon et al. 1984; see also Linz et 
al. 2008). If these objects were viewed down the outflow cav- 
ities, we would directly view the warm dust at the base of the 
cavity and perhaps from the accretion disk itself. This extinc- 
tion free view of the central regions would produce artificially 
steep laws from spherical models that cannot account for the 
geometry. 

The remaining objects have power indices that are shal- 
lower with power indices between p = 0.0 and p = 0.5. The 
best-studied of these objects is IRAS 20126. As discussed in 
Sect. 3.3.8 this object is a case where there is clear evidence 
that we are viewing the object close to edge-on. A dark lane ap- 
pears at wavelengths shorter than 20 pm. As a result, there will 
still be significant optical depth at 24.5 pm which suppresses 
the envelope emission. Instead the weaker, but more extended 
cavity wall emission becomes significant. This leads to an ap- 
parent need for shallower density profiles in our fits. 

We therefore conclude that in those MYSO cases where we 
view the objects at intermediate inclinations, the emission from 
the envelope dominates. Spherical models are then capable of 
reproducing simultaneously the 24.5 pm intensity profile and 
the SED. For these cases, radial density profiles with a power 
index of p = 1 .0 are preferred. 

5. Discussion 

We have presented resolved 24.5 pm images of a sample of 14 
MYSOs. In most cases, the MYSO is discrete, single (within 
< 2") and has a circularly symmetric profile on the sky. 
Simultaneous modelling of the spatial profile and the SED with 
simple spherical envelope models shows that those objects that 
can be adequately modelled in this way have radial density 
powerlaws n oc r~ p with exponent p — 1.0. 




0.0 0.5 1.0 1.5 2.0 2.5 3.0 

log(A) in micron 



Fig. 28. As Fig. 4, but for W3IRS5. The SED data is described in 
Sect. 3.4.3. The best fitting model has a p = 1.0 radial density profile. 



A number of studies have analysed the density structure 
of MYSO (and UCHu) envelopes by means of 1-D radiative 
transfer modelling. The approaches consist in reproducing the 
observed SED and the submm intensity profiles (Mueller et al. 
2002; Beuther et al. 2002; Williams et al. 2005) and molecular 
line emission (van der Tak et al. 2000; Hatchell & van der Tak 
2003). A few important differences exist between these studies 
and the present one. First, the COMICS observations provide 
information on scales about ten times smaller, i.e. 1000 AU. 
These scales could correspond to the transition region between 
the material reservoir (the envelope) and a putative accretion 
disk. Secondly, the above mentioned studies ease or ignore 
spectral constraints at wavelengths shorter than 100//m. In con- 
trast, our approach is to be more cautious with the IRAS mea- 
surements if they cannot be reproduced by any spherical model, 
but keep the mid-IR and submm measurements. This assump- 
tion derives from the likelihood of source confusion and back- 
ground contamination within the large IRAS beam, and our ex- 
perience with Spitzer 10pm fluxes which are often significantly 
lower than the interpolated IRAS fluxes, even after correction 
for detector non-linearity effects. In practice, this resulted in 
somewhat lower bolometric luminosities than usually adopted 
for our MYSO sample, primarily caused by the luminosity sen- 
sitive intensity profile. 
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Fig. 29. COMICS 24.5 fim mosaic of the Cep A region. Dark con- 
tours represent Spitzer 8/im emission. Centre of the coordinate sys- 
tem corresponds to the phase centre of the VLA 7 mm continuum map 
in Jimenez-Serra et al. (2007): a = 22 h 56 m 18.0 s , <S = 62°01'49.5" 
(indicated by a white cross). Contour levels of the 24.5 pm emis- 
sion are at 10%, 15%, 25%, 50%, and 80% of peak flux density 
(2.7 10 1 Jy arcsec~ 2 ). North is up, East is to the left. 

As far as the radial density distribution is concerned, we 
find a preference for p = 1 .0 models, as discussed in Sect. 4. 
The present analysis shows that radial profiles as steep as 
p = 2.0 are incompatible with the observations, and even less 
steep p = 1.5 cannot always be justified. Our values are consis- 
tent with the ones found by van der Tak et al. (2000) from the 
submm dust continuum and molecular line emission, albeit on 
much larger size scales. These authors find a range of 1.0 - 1.5 
for the powerlaw index. Williams et al. (2005) draw a similar 
conclusion for a large sample of candidate MYSOs. They find 
an average value of the powerlaw index of 1.3 ± 0.4. However, 
one should note that the models in the latter study require very 
high optical depths (generally optically thick at 100//m), which 
is a factor of a few to an order of magnitude larger than derived 
by us. It goes to show that the DUSTY models presented by 
Williams et al. would probably never fit the silicate feature nor 
the mid-IR part of the SED, and indeed they do not attempt to 
fit the mid-IR nor the far-IR data points. Mueller et al. (2002) 
claim a value for the power index "p" which is consistent with 
the previous two but slightly higher: 1 .8 ±0.4. Importantly, they 
argue that the van der Tak (2000) results would favour higher 
values for the power index if they had convolved their models 
with the actual telescope beam instead of a Gaussian profile. 
Finally, Beuther et al. (2002) claim 1.6 ± 0.5 from powerlaw 
fits to the inner 32" of their resolved 1.2 mm images. In sum- 
mary, the far-IR and submm studies generally prefer somewhat 
steeper density profiles than the ones derived from 24.5 /im 
images, recalling the caveats regarding the van der Tak et al. 
(2000) and Williams et al. (2005) results. This apparent incon- 



sistency can however be brought into agreement when allowing 
for the different spatial scales probed. 

On small angular scales, various studies in the (sub)mm 
find evidence for shallower density distributions or flat inten- 
sity distributions, which are suggested to be due to an unre- 
solved core or possibly a collection of cores. For example, 
Hatchell et al. (2000) find better model fits to their submm 
intensity profiles of a number of UCH n regions, when they in- 
clude an unresolved, high optical depth, central core. Beuther et 
al. (2002) describe the 1 .2 mm intensity profiles of a large sam- 
ple of massive cores with models that require an unresolved, 
inner, constant intensity distribution with radii between 2 000- 
60 000 AU. Van der Tak et al. (1999, 2000) conclude from their 
interferometric mm data (probing ~1000AU scales) that the 
detected emission is caused by a compact structure with an es- 
timated angular size of 0.3", and which is different from the 
larger scale spherical envelope structure. 

Compact substructures thus seem to be required in order to 
explain the (sub)mm observations. Flat intensity distributions 
are suggestive of an optically thick dust component, such as a 
dense shell or a disk (van der Tak et al. 2000), or a collection 
of subcores created by the ongoing fragmentation of the large- 
scale molecular core (Beuther et al. 2002). In the latter case, the 
low angular resolution single dish submm observations would 
find the integrated emission of the various subcores to be equiv- 
alent to a shallow or a flat density distribution. It can therefore 
be argued that the flattening of the radial density distribution 
seen in the submm close to the unresolved central source is con- 
sistent with the relatively shallow density distribution as traced 
by our high-resolution 24.5 yum images. 

The images are generally dominated by one single, resolved 
principal source, identified with the MYSO, which has a circu- 
lar appearance on the sky. Only in a small minority of cases 
do we find more than one condensation on scales < 2". This 
observation goes against the idea of fragmentation as the main 
cause of the flattening of the inner density profile. Rodon et al. 
(2008) present 0.35" resolution mm images of the near-IR dou- 
ble source W3 IRS5, in which the dominant component does 
not resolve into multiple mm sources. Given that this source 
is the most luminous source in our sample, any fragmentation 
into a cluster of protostars would be expected in particular in 
this source, but this is not observed. Instead, it seems that the 
24.5yum could indicate that it is the density distribution of the 
ambient core material itself that actually flattens out at scales 
of an arcsecond. An alternative explanation why this could be 
so is rotation. The density structure for a rotating and infalling 
envelope (TSC envelopes, Terebey et al. 1984) becomes on av- 
erage significantly shallower within the so-called centrifugal 
radius, i.e. where rotational motion dominates over infall. The 
evolution of the density distribution from relatively steep on 
large scales as seen by single-dish the (sub)mm observations 
to shallow on smaller scales at shorter wavelengths could be a 
manifestation of this effect. 

The inability of models with spherical and smooth den- 
sity distributions to reproduce the near-IR and mid-IR part of 
the SED was demonstrated again in the analysis presented in 
this paper. The development of more sophisticated models (e.g. 
Yorke & Sonnhalter 2002; Whitney et al. 2003; Indebetouw et 
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Table 5. Overview of the quality of the simultaneous model fits to 
the SED and the intensity profile. The second column gives the power 
index of the radial density distribution corresponding to the best fit- 
ting model. A colon is added to uncertain power index values. Model 
parameters are given in table 4. A tickmark indicates a reasonable cor- 
respondence between model and observations; an "x" means no cor- 
respondence, and indicates that not enough data or only data of 
insufficient quality is available. 



sive star forms can only be settled with such spatially resolved 
information. 



6. Conclusions 

We have presented a study of resolved 24.5 /im images of 
14 massive star forming regions. The images probe linear 
size scales of 1000 AU. Emission at 24.5 /im is dominated by 
th e MYSOs in these regions. They are discrete sources and 
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al. 2006) is driven in part to reproduce this emission. Radiative 
transfer models that incorporate the mentioned rotating TSC 
envelopes and implement a flared equatorial accretion disk and 
a bipolar outflow cavity are presented in Whitney et al. (2003). 
The geometrical features are inspired mainly by detailed obser- 
vations of low-mass class I YSOs (Whitney et al. 2003; see also 
Tobin et al. 2008) and may therefore be considered more real- 
istic than the simple spherical models we have adopted here. 

Recent work has applied the more sophisticated envelope 
models to MYSO SEDs and concluded that high-mass star for- 
mation proceeds similar to low-mass star formation (e.g. De 
Buizer, Osorio & Calvet 2005; Fazal et al. 2007). Such work 
is becoming popular especially after the publication of an ex- 
tensive grid of SEDs (Robitaille et al. 2006) that are based on 
the Whitney et al. models. The rationale behind the Robitaille 
et al. SED grid is that high-mass stars may form similarly to 
low-mass stars, or at least that the SEDs of MYSOs are deter- 
mined by similar geometrical structures as found in low-mass 
objects. Inferring 2D (or 3D) information from SED fits only is 
prone to be subjective and misleading. The quoted conclusion 
therefore that massive SF is similar to low-SF based on SED 
fits to such models is a circular argument. What the Whitney et 
al. models in particular show is that one should be careful inter- 
preting SEDs, without any independent additional data. Large 
scale spatial information of the dust emission in MYSOs can 
be retrieved from the single-dish (sub)mm studies previously 
mentioned, intermediate scales can be probed in the mid-IR 
as our study presented here, and with (sub)mm interferometry. 
Finally, scales down to 100 AU in the mid-IR can be reached 
with mid-IR interferometry using e.g. the VLTI (de Wit et al. 
2007; Linz et al 2008). Eventually the question of how a mas- 



velope. Various regions display extended, diffuse emission. 
This emission is associated with UCHii regions (e.g. Mon 
R2, AFGL437, AFGL2591) in which case the dust is heated 
by Lya within the ionized zone. Shock excited material (e.g. 
S 140) also seems to produce diffuse emission at 24.5 yum. 

Simple 1-D spherical model fits to the MYSO 24.5 yum spa- 
tial profile and SED show that radial density powerlaws of the 
form n - no (r/r su bi)~ p with a power p = 1 .0 are preferred. 
When there is evidence that we are viewing the MYSO either 
face-on down the outflow cavity or edge-on through a torus we 
find steeper or shallower density laws respectively . These den- 
sity laws are more likely to be due to the inadequacies of the 
spherical models in these cases than a true representation of a 
different density law, and we expect that p = 1 .0 also applies 
there to. 

We find that the spatial profile of the dust emission on 
scales of 1000 AU is shallower than that from larger 10 000 AU 
scales probed by (sub)mm dust emission. Inner flattenings seen 
in the submm are consistent with our results here. This flatten- 
ing is not likely to be due to fragmentation of the core, but 
due to the actual distribution of the emitting material. This 
is supported by the relatively small multiplicity of condensa- 
tions seen on sub-arcsecond scales. The continuous flattening 
from large to small scales could be the manifestation of rota- 
tion, but this requires further study. The application of more 
sophisticated, multi-dimensional models in relation to the data 
presented here and at even higher resolution using mid-IR in- 
terferometric observations will be the subject of a future paper. 
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Table 4. Overview of the modelling results. Graphical representation of the models are given in the figures to each subsection. Models #1 
are the best-fitting models. Models #2 and #3 are the best fitting models for radial density powerlaws that bracket the best-fitting density 
powerlaw of model #1 (except fot AFGL2591). Note that for steeper density profiles to fit, larger dust optical depths are required. Distances 
are taken from the listed references (Col. 2), and observed luminosities are given in Col. 3. Their respective references are given in Col. 4. The 
columns from Col. 6 and onwards are model parameters. Col. 7 is the power index of the radial density powerlaw of the model calculation, i.e. 
n = n (r/r sabl )~ p - Col. 8 is the dust type (see Sect. 3.2). Col. 11 is the ratio of the outer radius to the dust sublimation radius (r suW ). Col. 15 is 
the angular size on the sky of the inner dust free region. 
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